﻿using System.Drawing.Design;
using RayDen.Library.Core.Components;
namespace RayDen.Library.Core
{


    public partial class SobolSequence : BaseQMCSequence
    {
        private const int SOBOL_BITS = 32;
        private const int SOBOL_MAX_DIMENSIONS = 21201;

        private struct SobolDirectionNumbers
        {
            public int d, s, a;
            public int[] m;

            public SobolDirectionNumbers(int m_d, int m_s, int m_a, int[] m_M)
            {
                d = m_d;
                m = m_M;
                s = m_s;
                a = m_a;
            }
        }

        private uint[] directions;

        public SobolSequence(int? size)
        {
            RequestSamples(size??64);
        }
        void RequestSamples(int size)
        {
            directions = new uint[size * SOBOL_BITS];
            SobolGenerateDirectionVectors(directions, size);
        }

        private int pass, index;
        float rng0, rng1;
        public override float NextFloat()
        {
            index++;
            if (index == int.MaxValue)
                index = 0;
            if (index >1280000 && index%1280000 == 0)
            {
                pass++;
            }
            uint result = SobolDimension(this.pass, index);
            float r = result * (1f / 0xffffffffu);

            // Cranley-Patterson rotation to reduce visible regular patterns
            float shift = (index & 1)>0 ? rng0 : rng1;

            return r + shift - MathLab.Floor(r + shift); 
        }

        uint SobolDimension(int index, int dimension)
        {
            var offset = dimension * SOBOL_BITS;
            uint result = 0u;
            var i = index;
            for (uint j = 0;i != 0; i >>= 1, j++)
            {
                if ((i & 1)!=0)
                    result ^= directions[offset + j];
            }

            return result;
        }
        void SobolGenerateDirectionVectors(uint[] vectors, int dimensions)
        {
            Assert.IsTrue(dimensions <= SOBOL_MAX_DIMENSIONS);

            const int L = SOBOL_BITS;

            // First dimension is exception
            var v = vectors;

            for (var i = 0; i < L; i++)
                v[i] = (uint)(1 << (31 - i)); // all m's = 1

            for (var dim = 1; dim < dimensions; dim++)
            {
                SobolDirectionNumbers numbers = SOBOL_NUMBERS[dim - 1];
                var s = numbers.s;
                var a = numbers.a;
                var m = numbers.m;

                //v = vectors[dim * L];

                if (L <= s)
                {
                    for (var i = 0; i < L; i++)
                        v[i] = (uint)(m[i] << (31 - i));
                }
                else
                {
                    for (var i = 0; i < s; i++)
                        v[i] = (uint)(m[i] << (31 - i));

                    for (var i = s; i < L; i++)
                    {
                        v[i] = v[i - s] ^ (v[i - s] >> s);

                        for (var k = 1; k < s; k++)
                            v[i] ^= (uint)(((a >> (s - 1 - k)) & 1u) * v[i - k]);
                    }
                }
            }
        }
    }

}